#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <pthread.h>
#include "common.h"

//
//  global variables
//
int n, n_threads;
particle_t *particles;
FILE *fsave;
pthread_barrier_t barrier;

// The number of fast multipole method levels
int depth_level;
// This will be the array containing all of the boxes/simulation domains
box ***boxes;
// The root box
box world;

//
//  check that pthreads routine call was successful
//
#define P( condition ) {if( (condition) != 0 ) { printf( "\n FAILURE in %s, line %d\n", __FILE__, __LINE__ );exit( 1 );}}

//
//  This is where the action happens
//
void *thread_routine( void *pthread_id )
{
	int thread_id = *(int*)pthread_id;

	int particles_per_thread = (n + n_threads - 1) / n_threads;
	int first = min(  thread_id    * particles_per_thread, n );
	int last  = min( (thread_id + 1) * particles_per_thread, n );
	
	//--------------------------------------------------------------

	// This part will calculate which of the boxes this thread is responsible for
	// Each thread only calculates for a part of the simulation domain
	// This is done at level 2 where there are 16 boxes
	// at maximum this means we can have 16 threads
	int x = thread_id / 2;
	int y = thread_id % 2;
	
	int avg = 16 / n_threads;
	int mod = 16 % n_threads;
	
	int index;
	int next_index;
	
	if(mod == 0)
	{
		index = thread_id * avg;
		next_index = (thread_id + 1) * avg;
	}else if(thread_id < mod)
	{
		index = thread_id * (avg + 1);
		next_index = (thread_id + 1) * (avg + 1);
	} else
	{
		index = mod * (avg + 1) + (thread_id - mod) * avg;
		next_index = mod * (avg + 1) + (thread_id - mod + 1) * avg;
	}
	
	// Convert index into row and column numbers
	int l2x = index / 4;
	int l2y = index % 4;
	int l2xe = min((next_index - 1) / 4, 3);
	int l2ye = min((next_index - 1) % 4, 3);

	//----------------------------------------------------

	// box_list is a list each thread has, in the case of serial there is only 1 box_list
	// box_list contains the list of all the boxes, on each layer, that fmm must go through and calculate
	// This is done here so that fmm can just do calculations instead of figuring out in which
	// boxes to do the calculation
	vector < vector < pair < int, int > > > box_list;
	for(int i = 0; i <= depth_level; i++)
	{
		vector < pair<int, int > > dummy_v;
		box_list.push_back(dummy_v);
	}

	// This part wanders through the level 2 boxes this process got assigned to itself
	// It calls the function find_bottom_box(), this function goes to the bottom layer
	// children of this box and add them to this process's bottom_boxes list.
	int tcx = l2x;
	int tcy = l2y;
	for(int nn = next_index; index < nn; nn--)
	{
		if(tcy > 3)
		{
			tcy = 0;
			tcx++;
		}
		find_bottom_box(boxes[2][tcx][tcy], boxes, box_list, depth_level, tcx, tcy);
		tcy++;
	}
	
	//
	//  simulate a number of time steps
	//
	for( int step = 0; step < NSTEPS; step++ )
	{
		// Resets particle force
		for( int i = first; i < last; i++ )
			particles[i].ax = particles[i].ay = 0;

		pthread_barrier_wait( &barrier );
		
		// Calls Fast Multipole Method
		fmm(boxes, particles, depth_level, box_list);
		
		pthread_barrier_wait( &barrier );

		//  move particles
		for( int i = first; i < last; i++ ) 
			move( particles[i] );
		
		// If thread is 0 send root particles to level 1
		if(thread_id == 0)
			sort_elements(world, particles, boxes, 0, 0);

		pthread_barrier_wait( &barrier );

		// Send particles from level 1 to level 2
		if(thread_id < 4)
			sort_elements(boxes[1][x][y], particles, boxes, x, y);
		
		pthread_barrier_wait( &barrier );
		
		// Send the particles to the rest of the boxes
		assign_elements(boxes, particles, 2, depth_level, l2x, l2y, l2xe, l2ye);

		//  save if necessary
		if( thread_id == 0 && fsave && (step%SAVEFREQ) == 0 )
			save( fsave, n, particles );
	}

	return NULL;
}

//
//  benchmarking program
//
int main( int argc, char **argv )
{    
	//
	//  process command line
	//
	if( find_option( argc, argv, "-h" ) >= 0 )
	{
	printf( "Options:\n" );
	printf( "-h to see this help\n" );
	printf( "-n <int> to set the number of particles\n" );
	printf( "-p <int> to set the number of threads\n" );
	printf( "-o <filename> to specify the output file name\n" );
	return 0;
	}

	n = read_int( argc, argv, "-n", 1000 );
	n_threads = read_int( argc, argv, "-p", 4 );
	char *savename = read_string( argc, argv, "-o", NULL );

	//
	//  allocate resources
	//
	fsave = savename ? fopen( savename, "w" ) : NULL;

	particles = (particle_t*) malloc( n * sizeof(particle_t) );
	set_size( n );
	init_particles( n, particles );
    
	//-----------------------------------------------------

	// Creates the main box, this is box number 0, the root
	init_world(&world);
	
	// Adds the particles to the world box
	for(int i = 0; i < n; i++ )
		world.add_new(i);
	
	// Calculates what the depth level should be
	// The depth level is the number of layers our simulation will have
	depth_level = ceil(log10(n));
	printf("depth_level = %d \n", depth_level);
	
	// Each layer is a 2 dimensional array containing boxes
	// layer 0 has 1 box
	// layer 1 has 4 boxes
	// layer 2 has 8 boxes
	// etc...
	boxes = new box**[depth_level + 1];	// Creates the level 0 multi dimensional array

	// Initializes the 2 dimensional arrays, creates the boxes
	initializations(boxes, depth_level);
	// Links the root boxes 4 children to it and sends them their particles
	split(world, particles, boxes, 0, 0);
	// Links the 4 level 1 children with their own children
	// and finishes linking all the other boxes with their children
	// also, finishes distributing all the particles to each box accordingly to their coordinates
	build_tree(boxes, particles, 1, depth_level, 0, 0, 1, 1);

	//-----------------------------------------------------
	pthread_attr_t attr;
	P( pthread_attr_init( &attr ) );
	P( pthread_barrier_init( &barrier, NULL, n_threads ) );

	int *thread_ids = (int *) malloc( n_threads * sizeof( int ) );
	for( int i = 0; i < n_threads; i++ ) 
		thread_ids[i] = i;

	pthread_t *threads = (pthread_t *) malloc( n_threads * sizeof( pthread_t ) );

	//
	//  do the parallel work
	//
	double simulation_time = read_timer( );
	for( int i = 1; i < n_threads; i++ ) 
		P( pthread_create( &threads[i], &attr, thread_routine, &thread_ids[i] ) );

	thread_routine( &thread_ids[0] );

	for( int i = 1; i < n_threads; i++ ) 
		P( pthread_join( threads[i], NULL ) );
	simulation_time = read_timer( ) - simulation_time;

	printf( "n = %d, n_threads = %d, simulation time = %g seconds\n", n, n_threads, simulation_time );

	//
	//  release resources
	//
	P( pthread_barrier_destroy( &barrier ) );
	P( pthread_attr_destroy( &attr ) );
	free( thread_ids );
	free( threads );
	free( particles );
	if( fsave )
	fclose( fsave );

	return 0;
}
